Dynamics of Bose-Einstein Condensates in One-Dimensional Optical Lattices in the 

Presence of Transverse Resonances 
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The dynamics of Bose-Einstein condensates in the lowest energy band of a one- dimensional optical 
lattice is generally disturbed by the presence of transversally excited resonant states. We propose 
an effective one- dimensional theory which takes these resonant modes into account and derive vari- 
ational equations for large-scale dynamics. Several applications of the theory are discussed and a 
novel type of "triple soliton" is proposed, which consists of a superposition of a wavepacket at the 
upper band edge and two transversally excited wavepackets which are displaced in quasi-momentum 
space. 
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I. INTRODUCTION 

The phenomenon of Bose-Einstein condensation is a 
collective effect which relies on the bosonic nature of the 
particles alone (for reviews see, e.g., Refs. 0-0). Al- 
though an interaction between particles is not needed 
for the corresponding phase transition, its presence has a 
substantial influence on the properties of a Bose-Einstein 
condensate (BEC). In this context, solitons are of funda- 
mental interest since they represent states whose very 
existence relies on the interaction. 

For atomic BECs, bright solitons as well as dark soli- 
tons have been experimentally demonstrated for atoms 
with attractive 0, and repulsive interaction d [ljj , re- 
spectively. The present work is motivated by the recent 
observation of gap solitons in a 87 Rb BEC E2 Gap 
solitons are bright solitons for a BEC with repulsive in- 
teraction in an optical lattice and rely on the negative 
effective mass around the upper band edge of the peri- 
odic potential. To create a gap soliton it is necessary 
to control the motion of the initial wavepacket in quasi- 
momentum space This kind of physical situation 
has recentl y b een intensively studied, both theoretically 
El El E11I3 and experimentally El El Ell- 
in the present work we consider the influence of the 
transverse confining potential on the dynamics of a BEC 
in an one-dimensional optical lattice. We are particu- 
larly interested in the behaviour around the upper band 
edge of the lowest energy band. In this energy range the 
transverse confinement leads to the presence of transver- 
sally excited resonant states which significantly change 
the stability of the BEC [20L l2lj and alter its dynamics 
p2| . The resonances are important if the transverse exci- 
tation energy is small compared to the modulation depth 
of the optical lattice. 

Much of the recent research on BEC is concerned with 
an effectively one-dimensional situation. Generally this 
can be achieved if the transverse excitation energy is large 



compared to the interaction energy. This allows a sim- 
plified one-dimensional description of the dynamics by 
either projecting the collective wavefunction on the trans- 
verse ground state |2j, |6| or, more accurately, by making 
a Gaussian variational ansatz for the transverse shape of 
the wavefunction [2! |24| . While such an approach gives 
excellent agreement with a full three-dimensional theory 
in absence of transverse resonances (i.e., around the lower 
band edge in the case of a ID optical lattice [21|1. it is 
not suitable to describe a BEC around the upper band 
edge |51 |2lJ. In this paper we present a generalized 
one-dimensional theory by projecting the collective wave- 
function on a superposition of longitudinal wavepackets 
centered around the resonant states. In Sec. |n] we will 
review the preparation of a BEC at the upper band edge 
in order to motivate our particular approach. In Sec. IIIII 
simplified dynamical equations are derived and compared 
to previous approaches. In Sec. IIVI we further reduce 
these equations by making a variational ansatz for the 
wavefunction. In Sec. we will discuss several solutions 
of this system. 



II. DESCRIPTION OF THE PROBLEM 



Very recently, gap solitons have been experimentally 
observed in a BEC of Rubidium atoms Gap soli- 

tons correspond to a wavepacket of repulsively interact- 
ing atoms prepared at the upper band edge of the lowest 
band in an optical lattice. The process of creating a gap 
soliton is quite sophisticated since one has to move the 
BEC from the ground state, where it first is created, to 
the upper band edge of the optical lattice. For the pur- 
pose of this paper it can be summarized in the following 
way: first, a BEC is created in the ground state of a 3D 
harmonic trap Vtrap(x) = M 2 ivjz 2 /2+ M 2 ivj_(x 2 + y 2 )/2, 
where M is the atomic mass and u>i are the axial and 
transverse trapping frequencies. Then a one-dimensional 
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optical lattice of the form 

V opt {z,t) = V co${2k L z + 4>{t)) (1) 

along the z axis is switched on adiabatically and the ax- 
ial harmonic trap is switched off (wm = 0). Here, fcj, 
is the wavevector of the laser beam forming the optical 
lattice. At this time the lattice phase (j>{t) is zero. The 
BEC is thus prepared as a wavepacket around the lower 
band edge of the lowest band of the optical lattice. Fi- 
nally, a Bloch oscillation is employed (4>(t) varying with 
time) so that the wavepacket is slowly moving upwards 
in the energy band (so that excitations to higher bands 
can be neglected) and eventually reaches the upper band 
edge. This is an application of dispersion management 
for atomic matter waves which is described in more de- 
tail in Ref. |l2l | and is now of high experimental interest 

0EEI. 

To describe the dynamics of a BEC that is manip- 
ulated within the lowest energy band of the lattice, it 
would be desirable to have an effective dynamical equa- 
tion at hand which is one-dimensional and based on the 
effective-mass approximation, rather than including the 
full periodic and transverse trapping potentials. To de- 
rive such an equation we start from the Gross-Pitaevskii 
equation (GPE) for a BEC in a ID optical lattice and a 
transverse trapping potential, 

ihd t ip(-x, t) = (H\\ +H±- Maz) t) 

+ K |V(x,t)| 2 V(x,i) (2) 

with := (x,y) and 

H\\ = ^ + V opt (z,t) (3) 

H± = + Wx ± ). (4) 

Here tp is the collective atomic wavefunction which we 
assume to be normalized to one. The interaction param- 
eter is given by k := ^irh 2 a sc&t tN / M with a scat t being 
the atomic scattering length and N the number of atoms 
in the BEC. We have also included a homogeneous force 
term which corresponds to an acceleration a of the atoms. 
This term is closely related to the time variation of the 
lattice phase <p(t) and responsible for the generation of 
Bloch oscillations of the wavepacket. To avoid exciting 
atoms to higher bands of the optical lattice, the accelera- 
tion must be small enough so that adiabatic motion in the 
lowest band is possible. Throughout the paper we will 
assume that this is the case. We have omitted a longi- 
tudinal confining potential since our aim is to study the 
effects of the transverse dynamics rather than the per- 
turbation of the longitudinal lattice symmetry. A weak 
longitudinal potential could be included by introducing 
a slow variation of the lattice parameters |25|, however. 

Being nonlinear and inhomogeneous, Eq. © is impos- 
sible to solve analytically. Even the numerical simula- 
tion of it is time-consuming because of the necessity to 



resolve features on the scale of half the laser's wavelength 
(which is equal to the period of the lattice). In addition, 
it would be desirable to have a description which uses 
the (numerically verified) fact that the wavepacket stays 
localized in the energy band for a long time if the mod- 
ulation depth Vq is sufficiently small. We remark that, 
if Vq gets too large, a phase transition to a spatially lo- 
calized state which is smeared out over the lowest energy 
band takes place instead |2q . 

To derive such an analytical theory, we employ the ob- 
servation that a wavepacket, which is narrowly localized 
around a certain quasi-momentum go m the lowest en- 
ergy band, is very broad and varies slowly in position 
space. Let us assume for the moment that no transverse 
excitations are produced. Then one can make the ansatz 

ip{x,t) = B(z,t)ip qo (z)xo(x 1 _) (5) 

where ip qo is a quasiperiodic (Bloch) eigenfunction of 
H\\ with quasimomentum qq. The function xo denotes 
the transverse ground state of the trapping potential. 
The (dimensionless) function B(z, t) is an envelope which 
describes the large scale features of the wavefunction, 
whereas the small-scale features are included in ip qo . The 
basic idea of our approach is to average over the small 
spatial scales and to derive an effective equation for the 
large-scale behaviour of the wavefunction, i.e., for the 
envelope B(z, t). 

III. EFFECTIVE DYNAMICAL EQUATIONS 
FOR RESONANT MODES 

Before we can start to derive an effective equation, the 
ansatz © has to be generalized in two aspects. First, 
since our aim includes to describe the adiabatic Bloch 
oscillation from the lower to the upper band edge, we can- 
not assume that the quasimomentum is fixed, but have to 
admit a time-dependent qo(t). Secondly, we have to take 
into account transverse resonances which appear if qo(t) 
is in the vicinity of the upper band edge (see Ref. [2(J and 
Fig. Numerical investigations suggest that it is suffi- 
cient to include the two nearest resonances only, because 
all other resonances are negligibly populated. Therefore, 
the ansatz (JSJ needs to be modified to 

V>(x,i) = '%2B i (z,t)(p qt ( t )(z)x nt (x±), (6) 

i 

where, for our purposes, i runs from to 2, q{ denotes the 
quasimomentum around which each of the three modes 
is centered, and rij represents the transverse excitation 
number (no = 0, n\ = ri2 = 2 since, by symmetry, only 
even levels can be excited |2(j). -Bj is a slowly varying 
envelope function for each of the three modes. 

To derive an effective equation for the envelopes, we av- 
erage over the small spatial scale set by the lattice length 
L a = n/kL. Following standard methods, we introduce 
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an averaging function f av (z) which is slowly varying on 
the scale of L ai has a narrow support whose width cor- 
responds to the resolution of the effective equation, and 
which is normalized to one, J / av (z)dz — 1. The width 
of /av should be much smaller than the scale on which 
the envelopes Bi are varying. A function g(z) is then 
averaged by calculating ((g)) (z) := f dz'f av (z')g(z — z'). 
With this method, the envelopes can be extracted from 
the wavefunction by evaluating 



with the usual interaction mode integrals 



dz'f av (z') J d Xj_ Xn t (x±)Pq,( Z ~ Z ')tl)(xx, Z - z') 

^ S ni>nj J dz'f av (z)Bj(z - z')(p*.(z - z')ip qj (z - z) 

3 

Y / B 3 (z)S nt , n] I dz'U(z')<p* qi (z-z')y qj (z-z'), (7) 



where we have used that Bj is approximately constant 
over the support of / av and where the time dependence, 
for brevity, is dropped out. The integral in the last line 
can be evaluated as follows: for j = i the function |<p g J 2 
is periodic with period L a . Therefore, \ip qi \ 2 dz = 
L a /L since the Bloch functions are normalized (L is the 
quantization length). Since / av is roughly constant on 
the scale of L a , we find, by cutting the integral into bits 
of length L a , 

[ dz"f av (z-z")\<p qi (z")\ 2 a J2f^(z-mL a )^ 



1 

T 



(8) 



since the sum is just the discretized expression for a Ric- 
mannian integral over / av with dz — L a . For j ' i, 
consider first the case that the width Lf of / av is very 
large, Lt = L. Then the integral is simply the scalar 
product between the two modes and therefore zero un- 
less qj = qi. For sufficiently large Lf, the integral is still 
approximately zero if qi and qj are not too close to each 
other, since the product of the Bloch wavcfunctions then 
oscillates rapidly and averages to zero. Assuming that 
this is the case we find from Eq. (J7J 



«/ d 2 x xx * n y qi m z ) = 



Bi(z) 



(9) 



When we apply the same procedure (projecting onto the 
transverse modes and averaging over the longitudinal 
part) to the GPE and insert the ansatz ©, we are led to 

ihB t = fwxjnj + -)Bj +L^6 ni , nj ((y* Qi H\\ip qj Bj)) 



+K^B*B k B l l! j . kl I^. kl 

3,k,l 

+(q i -Ma)zB i , 



ij : kl 



dz<Pq i <Pq i <Pq k <Pqi 



1 ij;kl j a ^^XniXrijXnkXni 



(ii) 

(12) 



A dot denotes the derivative with respect to time. In 
the derivation of Eq. (|10[l we have exploited the fact that 
the averaging over the interaction integrals can be done 
in much the same way as for Eq. J7J: the averaged in- 
teraction integrals are again either periodic or rapidly 
oscillating and therefore do essentially acquire a factor 
1/L, which we have multiplied out in Eq. i|l(J|) . 

The last line in Eq. (|10(l deserves a comment. The 
homogeneous potential term —Maz simply survives the 
averaging procedure and is a direct consequence of the 
corresponding term in Eq. @ ■ The term proportional to 
qi arises from the time derivative on the left-hand side of 
Eq. (J2J which includes a term of the form qi(d qi ip qi )Bi. It 
is not hard to see that, provided the assumption that the 
wavepacket remains in the lowest energy band holds true, 
the derivative with respect to the quasi momentum can 
be approximated by d qi (p qi ~ iztp qi . The term is then of 
the same form as the homogeneous force and can be aver- 
aged in the same way. It is interesting to note that in the 
case of a simple Bloch oscillation caused by the homoge- 
neous force we have qi = Ma so that the linear potential 
is cancelled. This is nothing but a different description 
of the fact that a Bloch oscillation simply corresponds to 
a shift of a wavepacket in quasimomentum space, again 
under the condition that no higher bands are populated. 
This is the case for the main wavepacket in Fig. ^ for 
which the time dependence of its quasimomentum qo is 
simply a consequence of the induced Bloch oscillation. 
However, for the modes q\ and qi the time dependence 
of the quasimomenta is determined by a resonance con- 
dition and is not directly related to the Bloch oscillation. 
Hence, these two modes are subject to a renormalized 
homogeneous force. 

To perform the averaging over the longitudinal Hamil- 
tonian in Eq. i|l(J|) . we employ the well-known effective- 
mass method from solid state theory (see, e.g., Ref. 27]). 
Using that Bj(p q . is narrowly localized around quasi mo- 
mentum qi we can expand this expression in terms of 
Bloch wave functions ip qj +Aq > which are eigenfunctions 
of Hn with eigenvalues E(qj + Aq). The eigenvalue can 
be expanded to second order in Aq, resulting in 



H^Bj w J dAq (ip qj+Aq \ip qj Bj) 



(13) 



E(q 3 ) + v 3 Aq + j V qj +Aq 



(10) 



In this equation, we introduced two important physical 
parameters: the group velocity Vj := dE(q)/dq\ q=qj and 

Intro- 



the effective mass Mf := (d 2 E(q)/dq 2 \ q=q .)- 1 . 
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ducing the function 

Bj(z) := J dAqe iA *'{<p qj+ A 9 \<p gj B j ), (14) 

it is easy to see that the action of Hu can be expressed 
as 

(H^B^z) « J ^-(H^B^z') x 

J dAqe- lA ^'^ +Aq (z), (15) 
with the effective Hamiltonian 

H^=E( qj )- l hv 3 d z -^p. (16) 

This allows us to write the averaged Hamiltonian action 
appearing in Eq. l|lUfl in the form 

((^H^B,)) = J ^(H^B^z') J dz"f av (z") x 

J dAq e iA ^ z - z '- z "\*.{z - z") x 

U qj+ Aq(z - z"), (17) 

where u q are the periodic Bloch wavcfunctions, fq{z) — 
exjp(iqz)u q (z). Because of the averaging, we are inter- 
ested in distances z — z' much larger than L a . In this 
case, the phase factor in the integral over Aq varies much 
faster with Aq than the periodic Bloch function u qj+Aq . 
We therefore can replace the latter by u qj . The integral 
over Aq then becomes, on scales much larger than L a , 
the delta function 2i:d(z — z' — z") and we arrive at 

({flH^Bj)) = |dz"(i/ c ( g || B J )(z-z")/(z")x 
« (H^Bj)(z) J dz"f{z") x 

l f*q l { Z ") L PlA Z ") 

« K^H^Bj. (18) 

The last step in our derivation of effective equations 
for the envelope functions Bj is to show that Bi and 
Bi are, on average, equal. To do so, we first note that 
Bi = L((ip*.ip qi Bi)) since Bi is slowly varying. Inverting 
Ea. (|14fl . we can rewrite this as 

Bi(z) = L Jdz"f av (z") j dAq^z - z") x 

<p qt+ M(*-*") J ^ e ^ B &') ■ ( 19 ) 

It is then possible to repeat the argument given above for 
the action of Hk . Writing the quasiperiodic Bloch func- 
tions if q in terms of the periodic Bloch functions u q , we 



again find a rapidly oscillating exponential in Aq which 
results in a spatial delta function for large scales. Inte- 
grating this we find 

L<(|^ 4 | 2 B i )) = Lm\u q f)) 

* LBi({K\ 2 )) 

= Bi. (20) 

Using this identity we find for the effective equation de- 
scribing the large scale dynamics of the envelopes 

ihBi = tkJxim + ^Bi+H^Bi 

+ K B j B k B l I l;kl I ^M 

+ (q t - Ma)zBi . (21) 

For the case of a single wave packet centered around a 
fixed quasimomcntum, an equation similar to Eq. 121|) 
has also been derived using multiple-scale perturbation 
theory in the context of nonlinear optics |28l | and atom 
optics We have chosen a different approach 

since the inclusion of time-dependent quasi momenta is 
more obvious using the averaging method. In the follow- 
ing sections we will apply this equation to examine the 
conditions under which gap solitons can be formed and 
how they evolve in time. 

IV. DERIVATION OF VARIATIONAL 
EQUATIONS 

A major advantage of Eq. (|21[l. compared to the full 
GPE, is the simple form of the effective Hamiltonians 
H^g || . It describes interacting particles in a homogeneous 
external potential with different masses and velocities. 
This allows us to find a simplified analytical description 
and thus to gain more insight in the dynamics of a BEC 
in an optical lattice. Numerical simulations of the full 
GPE indicate that for each mode the wavepacket remains 
localized around qi for a long time if the optical lattice 
is not too deep. It is therefore reasonable to assume 
that the wavepackets can approximately be described as 
Gaussian wavepackets and to make a variational ansatz 
for them. Following the technique described in Refs. 0, 
132] ] , we first observe that Eq. (|2T)l can formally be derived 
from the Lagrangean 

i 

+ (Efa) + faxirii + i) - (Ma - qi)z\ \B,\ 2 

idz B; Bl B*d z B t ) + J^^Btf) 
+ f E B t B *i B kBill M ltyM- ( 22 ) 

i,j,k,l 



5 



A consistent variational ansatz for Gaussian envelopes is 
achieved by setting 



Bi(z,t) 



Mt) 



■. exp 



(z - Zl {t)f 



v / 7 r 1 /2 Wl (t) V 2 Wi {ty 
+ifc(t)z + i~ fl (t)z 2 ^ 



i4>i(t) 
(23) 



This describes a wavepackct of width Wi and amplitude 
Ai (having dimensions of length 1 / 2 so that Bi is dimcn- 
sionless). It is spatially localized around Zj and has 
an instantaneous energy of Ti4>%- Its mean velocity and 
its variance are given by (vi) — {fii + 2^iWi)/Mf s and 
A Vl = ^/2~f 2 w 2 + l/(2wf)/Mf. 

Inserting this ansatz for the envelopes in the La- 
grangean and extremizing the corresponding action in- 
tegral, we derived a set of 18 equations which describe 
the evolution of the three Gaussian wavepackets involved. 
This task, as well as the algebraic manipulations follow- 
ing below, are rather tedious and therefore have been 
completed using Mathematica ■ Since the variational 
equations are somewhat lengthy we exploited the special 
features of our system to reduce its complexity. To do 
so, we restrict our considerations to the case when the 
wavepackets are already at the upper band edge so that 
the quasi momenta are time-independent and given by 
qo = hk^ and q\ = hki, — Sq as well as 52 = hki, + Sq, 
where fcz, is the wavenumber of the optical lattice which 
appears in the optical potential Q. <5g is identical to 
12 — Qo- It can be derived from the resonance condition 
that the three energies + haj±(ni + |) for i = 0, 1, 2 
are equal. Setting this energy to zero we can also omit 
the corresponding terms in the Lagrangean. Because the 
wavepackets are already at the upper band edge we will 
also not need the homogeneous force to induce Bloch os- 
cillations, i.e., we set a = qi = 0. 

The special values of the quasi momenta imply that 
most of the interaction integrals lfj. kt are zero or have 
an identical value. This can be seen by expanding the 
Bloch wavefunctions in terms of momentum eigenstates, 
(Pq(z) = J2i c /(<?) exp(iz(q + 2lhk]j)). By Fourier trans- 
forming the stationary Schrodinger equation H»tp g — 
E(q)(p q , one finds the following equation for the expan- 
sion coefficients ci(q), 



relative to each other. In the limit of an infinite op- 
tical lattice, the interaction integral /J' . fei will therefore 
vanish if these phase factors do not exactly cancel each 
other. For instance, iJo-oi = because its integrand is 
proportional to exp(iSqx), but 7q . 12 7^ 0. This, in com- 
bination with Eq. (|25|) . ensures that all interaction inte- 
grals, except 4' ;oo and ^ii 5 ii = 4a ; 22 = ^i2 ; i2 as wel1 as 
4>i-oi = ^02 02 = 4>o-i2' d° van ish (in addition, the sym- 
metries L 



ij:kl 



Ihkl a '"l %lk 



Hj-ki nave to be taken 



into account). Thus, there are only three independent 
interaction parameters which we will denote by 



Kq 
K 1 



K 
K 



rll T± 

- , 00;00 J 00;00 > 



rll T± 

i ll;ll i ll;ll i 



rii r- 1 - 

i—t- J 01;01- , 01;01 
y'TTfl 



(26) 



Even with these simplifications the resulting equations 
are still very lengthy, but they admit the analysis of sym- 
metric solutions. By symmetry, we have vq = and 
t>2 = —V\ for the group velocities of the wavepacket, and 
Mf = Mf for the effective masses. Under these con- 
ditions one can show that zo = and (3q = are so- 
lutions of the variational equations. This result is intu- 
itively clear and just means that the central wavepacket 
remains at the upper band edge with mean position and 
velocity zero. In addition, symmetry implies that the 
two transversally excited wavepackets should evolve in 
an identical way, but with opposite mean velocities (be- 
cause their group velocities differ by a sign). We there- 
fore can set A2 — A\, 72 =71, 4>2 — 4>ii W2 = t«i, and 
02 = Z2 = —Zi, which reduces the number of in- 

dependent variational parameters to ten (four for go and 
six for q\ ) . In addition, the conservation of the number 
of atoms implies the constraint A^ + 2A\ — L, so that 
we are effectively left with nine independent parameters 
. The resulting variational equations are given by 



An = 



^ SW An A\ 2 Kp\ 
WqWi 



(27) 



E(q)ci(q) 



{q + 2lhk L f 



2M 



C M) + y(cz+i(g) + Q- 



(24) 

This equation shows that the expansion coefficients are 
real and that, if ci(hki, — Sq) is a solution, then so is 
ci{hkL + Sq) = c-i-i(hkL — Sq). Thus, we have the rela- 
tion 



(25) 



It is well known, and can be seen from the above ex- 
pansion, that Bloch wavefunctions are periodic up to a 
phase factor exp(iqx). Therefore, the three wavefunc- 
tions <p qi are oscillating with a phase factor exp(±i5qx) 



Zi = Vi + 



h/3i 2/U171 e *"i S" (1) An 2 zi K01 



Mf 



Mf 



Wq Wl 



> (28) 



2 h wq 70 e 

WO = h 



W KOI 



Mf 



Wo Wl 



(29) 



The amplitudes Ai are normalized to L because the full wave- 
function Biip qi should be normalized to one and ip Qi carries a 
factor of 1/V~L because of its normalization. 



6 



-i = ^e^+ (30) ft = -]^I- 2 ^i-^^+ (35) 



V (-S (3) +S (1 W-2S«zi 2 ) K(u e ^ V2A! 2 Z! ( Wl 2 +4z! 2 ) Kl V 



8e~^ Ai 2 zi 3 kqi 



wo wi 



W q 5 wi ^ In these eq uations we have introduced the notation w := 

y/wo + w\ and 



4e~^ A x 2 (w 2 -2zi 2 ) n 01 



2 



G 



-2z(0 o — 0i) 

C.C. 



( n/2 

2i(7o - 7i ) , 



* - srv - + + (32) c := , ; , f — ^ + " ■ (36) 



e ™? Ai 2 (wi 2 - 4zi 2 ) ki 
V2wi 5 

2n m e~^ Ao 2 (w 2 - 2zi 2 ) 



fe + ^- 2l (7o-7i))' 



The functions depend on fa, 7,, and u»i and do van- 
ish for fa — fa = 0. 



V. SPECIAL SOLUTIONS OF THE 



-3, 9 , m , „ „ N , VARIATIONAL EQUATIONS 

Koie (-C( 3 )+CW( W i 2 -2z! 2 )) 

1,70 Wl Initially empty transverse excited modes: A surprising 

consequence of the variational equations can be seen im- 
mediately: it follows from Eq. 127|) that, when all atoms 
K o , /„„are in the central wavepacket (Ai — 0), the amplitude Aq 



+ " " u + (33] 



2My S wo 2 4v / 2u'o and therefore also ^! will not change in time. Thus, the 

z 2 transversally excited wavepackets would never be pop- 

e Ai 2 (— C^ 3 ^ +3C^^wo 2 ) «oi ulated. This prediction is a direct consequence of the 

2 wo 3 wi ^ assumption lJo;Oi ( = ^00:02) = an d m striking con- 

z z tradiction to the numerical results of Ref. |20|. This dif- 

2e 'a 1 A\ (3wo 4 + 2ioi 4 + wo 2 (5wi 2 — 2zi 2 )) ^oiference can be resolved when one recalls the conditions 

• • II 

w 5 under which our analytical theory is valid. Jqo'OI = 

is exactly fullfillcd only in the limit of an infinite opti- 
cal lattice. In a finite lattice the fact (discussed above) 
ft % Zl 2 that the integrand is oscillating with a phase factor of 

2 Mf ff w\ 2 ~ 2 Af° ff u>i 4 Vl ^ exp(±i5qz) only leads to oscillations of Iq . x, 80 tnat ^ 

2 2 z 2 is zero on average only. Since our wavepackets have a 

^ft _|_ K ?Ai {5(1 + 2 e _2 "?) w 4 + finite width in quasimomentum space, there will be a fi- 

2Mf 4V2 Wl 5 V v ; 

nite excitation probability even when A\ = initially. In 
„ *? „ z ? . addition, our theory assumes that the three wavepackets 

( — 2" \ 2 2 — 3" 4 1 ... • 

— 1 + 2 e *"i J W\ z\ + 16 e "121 I + are not overlapping in quasimomentum space, since only 

2 under this condition the averaging method can yield rea- 

e~"7 A 2 k sonable results. In practice, this is not exactly fulfilled 

2zi 2 — wi 2 ) + and will lead to corrections to the prediction of the av- 

Wl eraged equations. However, the time scale for transverse 

C 1 (3wi 4 + 4zi 4 )) + excitation out of a central wavepacket is quite large (typ- 

_ *? 2 ically about 70 ms [20^ so that the averaged equations 

e Ap kqi . ^ 4 ^ g Wl 4 + 4 Zl 4 + should provide a valid description for shorter times. In 

w 5 fact, the present considerations may provide another rea- 

wq 2 (5 wi 2 + 2 zi 2 )) , son for the long time scales for transverse excitations. In 
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addition, during the preparation of the wavepacket at the 
upper band edge through Bloch oscillations, the transver- 
sally excited modes are populated. Therefore, an initial 
condition with A\(0) ^ is realistic when we describe a 
system that already is prepared at the upper band edge. 

On the other hand, when using the initial condition 
Ai(0) = we are left with a theory for the central 
wavepacket only, since there are never any transversally 
excited atoms to interact with. In this case our descrip- 
tion reduces to the case considered in Ref. (but with 
a negative effective mass) so that one can transfer most 
of the results to our case. We therefore will not discuss 
it further. 

Case of three initial gap solitons: Another case of in- 
terest is the case when all three wavepackets are initially 
forming independent gap solitons. That is, in the absence 
of mutual interactions each of the three envelopes corre- 
sponds to a stationary solution of the variational equa- 
tions with self-interaction. We can find these solutions by 
setting «oi = and removing the terms proportional to 
m exp(—2z 2 /wf), which describe the interaction between 
wavepacket q\ and qi (see above). It is easy to see that 
in this case the soliton solution is given by 7; = /5j = 
and zf o1 = v\t as well as 



II! 



sol 



V2h 



A s ° mM, 



off 



C 1 = M0) + 



3Af 



4a/2 



(37) 
(38) 



The question remains whether this solution is stable 
against the presence of the mutual interactions of the 
three gap solitions. To answer it, we have made a sta- 
bility analysis by linearizing the variational equations in 
the deviations from the soliton solution We set 

Wi = wf° l + eSwi (and similarly for the other dynamical 
variables) and consider all equations only to first order 
in e, whereby the mutual interaction terms are treated as 
of first order in e. This is justified since these terms all 
include a factor which exponentially decays in time and 
thus have limited influence. Such a factor arises because 
the three wavepackets all have different group velocities 
and thus separate after a short time, the exponential be- 
ing a consequence of the overlap between the Gaussian 
wavepackets. The resulting linearized equations are given 
by 



5A 



t^y 

V „„sol / 



2 sin(Ac/>) A s ° l Af 1 ' 



«oi 



Sw 



-V2A$> } ' 



5 jo ko w s ' 



sol^ 



(39) 



(40) 



A ° l2 6w oKo , 5A o1 5A k -(^2r) a . ... 

+ e "i U {t), 



4V2< o1 



,80l 



(42) 



A\° l2 {8f3 l +2t Vl <E 7l ) Kl wl° l -(^kT 

' e 1 /ziW> 



V2 



(43) 



5wi = -V2Af 2 5 7 i K% wf 2 + e V "'? o1 ; f Wl (t), (44) 



KiVit (Af^Sw! + 2Af 1 6A 1 

V2 w s ° l4 

-(^iot) 2 
-2u 1 57 1 + e fait), 



(45) 



. Af l Ki (A? i Sw 1 + 2 5A 1 w? 1 ) 
6li = ^ .4 + e 1 /tiW. 



2v / 2<' 
A s ° l2 5w! ki (2t 2 Vl 2 + w\° 



(46) 



AV2wf 



sol' 1 



+ 



Af 1 5 Ax ki (-2 1 2 Vl 2 + 5 tof 



2-y/2u>i 



sol- 



_(±^y 



(47) 



The functions f a (t) depend on the soliton solution pa- 
rameters and increase at most polynomially (degree less 
than 4) in time. They represent inhomogenities, simi- 
larly to the right-hand side of Eq. (|39|l . Because of the 
exponentially decaying factor, these terms are only im- 
portant for times t < wf /v±. Therefore, to analyze the 
stability of the soliton solution, it is sufficient to solve the 
homogeneous linearized equations for a general set of ini- 
tial conditions, since for large enough times this correctly 
describes the general solution. 

To reduce the length of the linearized equations we 
have made an additional approximation. Our numerical 
simulations of the full GPE indicate that, after the BEC 
has been transferred to the upper band edge, the num- 
ber of atoms in the go wavepacket is considerably larger 
than in the other two modes 2 . Since A 2 = LNs, where 



yjsol 2 fi w ^ yl so ' 8 A Kq ~ ( " " 1 ) 2 2 r ^' le var ' a t' ona l equations presented in this work would predict 

(5^ = — _ ® _| 2 2 — --|-e ffa{t)j (41) that all atoms remain in the transverse ground state since the 

2 \/2l/)n ' \[2 Wn°^ excited modes are initially (almost) empty. 
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A^.j is the initial number of atoms in each mode, one can 
see that Af l < A s Q ° l and therefore wf l > w s ° l . Assum- 
ing that this is the case, we here present the linearized 
equations only to second order in the ratio A s ° 1 /Aq° 1 . 

The general solution of the homogeneous linearized 
equations is not hard to find. One immediately sees that 
SAo and therefore, because of atom number conservation, 
also 5Ai are constant in time. Swq and £70 are then cou- 
pled to each other only so that Eqs. (pljl and (j^T)l are 
easily solved. Swo and £70 then generally show a purely 
oscillating behaviour. This solution can then be inserted 
into Eq. I|42l) for the homogeneous phase factor. The lat- 
ter then grows in time, in addition to some oscillating 
factors, proportional to 3tKoA s ol 5A (0)/(V2w s o1 ). When 
this expression is compared to the evolution of the soli- 
tion phase factor l|38[l it becomes obvious that this linear 
increase in 5<f>o just corresponds to a small deviation, pro- 
portional to 5Ao(0)/Aq o1 , from the unperturbed energy 
of the soliton. We therefore have shown that the cen- 
tral soliton around quasi momentum qg is stable against 
the interaction with the other two wavepackets since its 
stability does also not depend on the evolution of the 
deviations in these wavepackets. 

The situation is quite different for the transversally ex- 
cited modes. Repeating the steps leading to the solution 
for the central wavepacket, one can see that the solution 
for 5 Pi is given by 



6/3i(t) = 5/3 1 (0)-2« 1 tcos(n 1 t)«y7 1 (0) 

Vxtsm{Slit) / Mi(0) 5t0i(O) 



(48) 



Af 1 



itr 



sol 



with fii := (4/3)d0i ol /di. This growing oscillatory 
behaviour clearly indicates instability against any ini- 
tial deviations 8wi(Q), <L4i(0), 5-fi(0), which unavoidably 
are introduced by the interaction between the three 
wavepackets. 

It is worth to examine the origin of this instability more 
closely. Our arguments are based on the fact that the 
two transversally excited wavepackets move away from 
the central wavepacket. This happens because we have 
set fli = /? 2 = for the excited wavepackets, so that 
they propagate with the group velocity ±ui. Hence, af- 
ter some time the wavepackets are separated, so that the 
mutual interaction disappears and cannot cause insta- 
bility anymore. However, setting = in absence of 
mutual interactions creates another source of instability: 
even in a strictly one-dimensional situation, a gap (or 
bright) soliton with non-vanishing group velocity is only 
stable 3 if the phase factor exactly matches the group 
velocity, /3, = —Mf s vi/h . Therefore, the instability of 
the transversally excited wavepackets is the same as that 
of an isolated gap soliton with the wrong phase factor. 



3 For other values of (3i a wavepacket characterized by Eqns. 1371 
and 1381 is still a stationary solution of the nonlinear Schrodinger 
equation, but it is not stable. 



The only possibility to avoid this kind of instability 
is to choose the appropriate phase factors P2 = —0i = 
-Mfvi/h. As a consequence, the excited wavepackets 
would remain at their original position so that the mutual 
interaction would not decrease. Since the latter is a res- 
onant coupling between the three wavepackets a general 
superposition of three gap solitons would not correspond 
to a stationary solution anymore. In the next section we 
will demonstrate that for a particular choice of parame- 
ters this problem can be overcome. 



VI. TRIPLE SOLITONS 

A particularly interesting situation appears when one 
tries to construct stationary wavepackets which remain 
spatially localized around Z\ = 0. As is evident from 
Eq. (J2H1), this is only possible for j3 x = -Mfvi/h. In- 
terestingly, this condition also guarantees the validity 
of 0i = in Eq. (|35[l . so that this requirement is self- 
consistent. The remaining equations will only lead to a 
stationary solution if the populations of the three modes 
are constant, i.e., if Aq = 0. Apart from the trivial so- 
lutions Aq = or Ai = this can be achieved by the 
condition S^ 1 ' — 0. A natural solution to this condition 
is (j>i — <po = and 70 = 71 = 0, whereby the latter as- 
sumption also ensures that the widths of the wavepackets 
remain constant. A necessary condition for this to hold 
are the equations 



7o = 7i = 



(49) 



Using Eqns. (|31|) - l|34[l this leads to algebraic condi- 
tions on the widths and populations of the three modes. 
The simplest way to solve these algebraic conditions is 
to, first, fix the ratio between the widths according to 
W2 = T]Wi, where 77 is some positive number. In addition, 
we write Ki = NRi/L so that re, is independent of the 
total number of atoms N and remains finite when the 
quantization length L goes to infinity. For these settings 
we derived solutions of the algebraic conditions which de- 
termine N, wi, and the population distribution among 
the modes as a function of 77, Ri, Mf , and Vi. A particu- 
larly nice example is the case when all three wavepackets 
have equal width, wi = wq- The solution then becomes 
very compact and is given by 



N = 

Wl = 



= L- 



ZMfRi 



6M, off 



K oi 



2M ofi (K - 3/soi) + 3M 1 eft («i - 2R m ) 
(f -R m )+Mf(Ri-2R i)) 



2i; 1 (2Af cff (f 



(6/4 - R Ri)^j2>Mf{Mf 

h 



3(M ofi - Mf 1 ) 



-Mfvi 



-2Af cff 



(50) 



(51) 



(52) 



with W2 = w 1 — wq and := KiL/N being independent 
from the number of atoms and the quantization length. 
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The state characterized by Eqs. Q50 [l -Q52 [l . which we 
will refer to as "triple soliton" , represents a special co- 
herent superposition of a wavepacket in the transverse 
ground state at the upper band edge of the optical lattice, 
and two wavepackets around the transverse resonances. 
The special choice 1)50 (1 -1)52 ^ for the parameters ensures 
that the mutual and self-interaction of the wavepackets 
exactly cancel the dispersion of each wavepacket due to 
its negative effective mass. It also guarantees, within the 
approximation that only two resonances are taken into 
account, that the triple soliton does not spread in the 
transverse direction. It therefore can be seen as a gen- 
eralization of the gap soliton which is unstable against 
transverse decay. It differs from the case of a superpo- 
sition of three gap-soliton wavepackets discussed above 
in that the mutual interaction between the wavepackets 
destroys the latter. This is because the stability criterion 
()37[) and ()38[l takes only into account the self-interaction 
of each of the three wavepackets. For the triple soliton 
the mutual interaction is included as well. 

A very interesting feature of the triple soliton is that 
the width of the soliton does not depend in any way on 
the interaction parameters of the system. It is solely de- 
termined by the structure of the lowest energy band of 
the optical lattice and in particular is proportional to the 
de Broglie wavelength of a particle of mass — Mf s mov- 
ing with the velocity v\ . The number N of atoms in the 
soliton depends on the interaction parameters, but it van- 
ishes if the group velocity v\ of the transverse resonances 
goes to zero, i.e., if the resonances are close to the band 
edge. The population of the three modes depends on the 
interaction and leads to consistency requirements: Since 
Aq can only take values between and L we find that 
the soliton can only exist if the effective masses fulfill the 
inequality 



3k i 



cir 



< 



Ml 



< 



2k, 



01 



(53) 



To see if this condition can be fullfillcd, we have nu- 
merically calculated the band structure for a BEC in 



a periodic potential of the form Vb cos(2fc£z), where 
ki = 27r/A^ is the laser's wavenumber and Vq the depth 
of the optical lattice, which we will give in units of the re- 
coil energy En = (l/2)Mv^ with the recoil velocity vr — 
hk L /M. We consider 87 Rb atoms (M = 1.45 x 1(T 25 kg, 
flscatt = 4.9nm) in an optical lattice driven by a laser 
close to the D2 line (\l = 780nm) and a 2-dimensional 
transverse harmonic trap of strength ui — 534s -1 . The 
effective mass, the group velocity, and the interaction 
parameters as a function of Vq arc shown in Fig. |2 a) 
and Fig. [3J respectively. As can be seen from Fig. [21 b) 
condition l|53|) can be fulfilled in this parameter regime, 
which also lies well within the range of current experi- 
ments 0, El ■ In Fig. 01 we have plotted the width as 
well as the number of atoms and population distribution 
for the novel kind of soliton. For the case wq = w\ under 
consideration, the population in the transversally excited 
modes is larger than that of the central wavepacket. 
VII. CONCLUSION 



Using an averaging method we have derived effective 
field equations which describe the large-scale behaviour 
of a transversally confined BEC in a one-dimensional op- 
tical lattice. Due to the existence of transversally excited 
modes resonant to wavepackets in the transverse ground 
state, these equations have the structure of coupled one- 
dimensional particles with different effective masses and 
dynamical interaction parameters. We have made a 
Gaussian ansatz for the envelopes of a wavepacket pre- 
pared at the upper band edge and the two nearest reso- 
nances in quasi-momcntum space. Variational equations 
for this ansatz are derived and several solutions are dis- 
cussed, including a novel kind of "triple" soliton. 
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FIG. 1: Scheme of the collective wavepacket's motion through 
the lowest energy band. The dotted lines represent the spec- 
trum of noninteracting atoms in a ID optical lattice and a 
transverse harmonic trap. The lowest of these lines corre- 
sponds to the lowest energy band of the lattice for atoms in 
the transverse ground state. The two upper copies of it are 
transversally excited atoms in the same band. The BEC is ini- 
tially prepared as a wavepacket around the lower band edge 
(lower left corner) and is adiabatically moved to the upper 
band edge with quasimomentum qo (dashed arrow). Around 
the upper band edge transversally excited resonances occur 
at quasimomenta q\ and 52- 



a) 

IM^'I/M, Vl /v R 




FIG. 2: a) Group velocity and absolute value of effective 
masses as a function of the optical lattice depth Vq. b) Fullfill- 
ment of condition 15311 as a function of Vo. The solid line rep- 
resents M° S /Mq S , the dashed lines are the upper and lower 
bound in the inequality 15311 . For Vo > OAEr the condition 
is fullfilled. 
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FIG. 3: Interaction parameters as a function of lattice depth 
Vo in units of the recoil velocity vr = ^/2Er/M. Solid line: 
K0j dashed line: Ri, dot-dashed line: Rqi- 
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FIG. 4: a) Width of the soliton wavepackets as a function of 
Vo. b) Total number of atoms N (solid line) in the soliton, 
and number of atoms Ni = A?N/L in mode i = (dashed 
line) and i = 1 (dotted line), respectively. 



